R Notebook

# clean up R environment
rm(list = ls())
library(dplyr)
library(tidyr)
library(stringr)
library(mosaic)
library(plyr)
library(ggplot2)
library(rvest)
library(reshape2)
library(readxl)
library(rpart)
library(partykit)
library(mclust)
library(tibble)
library(readr)
library(viridis)
library(RColorBrewer)
library(maps)
library(maptools)
library(mapdata)
library(ggthemes)
library(cluster)
library(fpca)
library(doParallel)

#using multiple cores
myCluster <- makeCluster(detectCores() - 1, type = "FORK")
# options
# display as decimal, not scientific notation
options(scipen = 999)
# data for 2013 - 2016
county_poverty <- read.csv(file = 'county_poverty_hist.csv', header = TRUE, sep = ',')
county_poverty$geo_sumlevel <- gsub('[[:punct:]]\\s', '', str_extract(county_poverty$geo_name, '\\,\\s[A-Z]{2}'))
county_poverty$geo_name <- gsub('\\,\\s[A-Z]{2}', '', county_poverty$geo_name)
county_poverty <- county_poverty[-3]
# naming scheme:
# pYEAR = number of people in poverty, for the given year
# prYEAR = poverty rate, for the given year
# popYEAR = total population, for the given year
names(county_poverty) <- c('county', 'state', 
                           'p2013', 'p2014', 'p2015', 
                           'pr2013', 'pr2014', 'pr2015', 
                           'pr2016', 'p2016', 'pop2013', 
                           'pop2014', 'pop2015', 'pop2016')
county_poverty <- county_poverty %>% 
  na.omit() %>%
  select(sort(current_vars()))
# remove Puerto Rico
county_poverty <- county_poverty %>% 
  filter(state != 'PR')
county_poverty <- county_poverty %>% 
  select(sort(current_vars()))

# economic/unemployment data
# https://apps.bea.gov/regional/
ec <- read.csv('CAINC4__ALL_STATES_1969_2017.csv', header = TRUE, sep = ',')
ec <- ec[, c('GeoFIPS', 'GeoName', 'Description', 'X2013', 'X2014', 'X2015', 'X2016')]
ec$Description <- gsub('[[:punct:]]', '', ec$Description) %>% 
  trimws()
# keep 'Per capita personal income (dollars) 4/', 'Total employment', 'Farm proprietors' income', 'Nonfarm proprietors' income'
# clean data
ec <- ec %>% 
  filter(ec$Description == 'Per capita personal income dollars 4' | 
                      ec$Description == 'Total employment' | 
                      ec$Description == 'Farm proprietors income' | 
                      ec$Description == 'Nonfarm proprietors income')
ec$Description[ec$Description == 'Per capita personal income dollars 4'] <- 'income'
ec$Description[ec$Description == 'Total employment'] <- 'empl_count'
ec$Description[ec$Description == 'Farm proprietors income'] <- 'farm_income'
ec$Description[ec$Description == 'Nonfarm proprietors income'] <- 'nonfarm_income'
ec <- dcast(melt(ec, 
                 id.vars = c('GeoFIPS', 'GeoName', 'Description'), 
                 measure.vars = c('X2013', 'X2014', 'X2015', 'X2016')), 
            GeoFIPS + GeoName + variable ~ Description) %>% 
  mutate(variable = gsub('X', '', variable))

attributes are not identical across measure variables; they will be dropped

ec[, 1:2] <- lapply(ec[, 1:2], as.character)
ec[, c(1, 3:7)] <- lapply(ec[, c(1, 3:7)], as.numeric)

NAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercion

ec$p_farm_income <- ec$farm_income / ec$nonfarm_income
ec[is.na(ec)] <- 0
# combine economic data with land area
area <- read_excel('LND01.xls')[,c(2,4)]
names(area) <- c('GeoFIPS', 'area')
area <- lapply(area, as.numeric)
ec <- merge(ec, area, by = 'GeoFIPS')
# state data
## COMBINE WITH 'STATE_DATA' DATA FRAME
state_ec <- ec[(ec$GeoFIPS %% 1000 == 0), ]
state_ec <- state_ec %>% 
  filter(GeoFIPS < 90000 & GeoFIPS != 0)
state_ec <- subset(state_ec, select = -c(GeoFIPS))
names(state_ec) <- c('state', 'year', 
                     'empl_count', 'farm_income', 
                     'income', 'nonfarm_income', 
                     'p_farm_income', 'area')

# create tall data set
cp1 <- gather(county_poverty[c(1:5, 14)], year, poverty_count, p2013:p2016, factor_key = TRUE) %>% 
  mutate(year = str_extract(year, '\\d{4}'))
cp2 <- gather(county_poverty[c(1, 6:9, 14)], year, total_count, pop2013:pop2016, factor_key = TRUE) %>% 
  mutate(year = str_extract(year, '\\d{4}'))
cp <- merge(cp1, cp2, by = c('county', 'state', 'year')) %>% 
  mutate(poverty_rate = poverty_count / total_count,
         year = year %>% as.numeric())
cp$county <- gsub('\\sCounty', '', cp$county)
cp$county <- gsub('\\sParish', '', cp$county)
cp$county <- gsub('\\scity', '\\sCity', cp$county)
cp$county <- gsub('[[:punct:]]', '', cp$county)
cp$county <- cp$county %>% 
  trimws()
# county data
county_ec <- ec[!(ec$GeoFIPS %% 1000 == 0), ]
county_ec <- subset(county_ec, select = -c(GeoFIPS))
county_ec$state <- gsub('[[:punct:]]\\s', '', str_extract(county_ec$GeoName, '\\,\\s[A-Z]{2}'))
county_ec$GeoName <- gsub('\\,\\s[A-Z]{2}', '', county_ec$GeoName)
names(county_ec) <- c('county', 'year', 'empl_count', 
                      'farm_income', 'income', 'nonfarm_income', 
                      'p_farm_income', 'area', 'state')
county_ec$county <- gsub('[[:punct:]]', '', county_ec$county)
county_ec$county <- gsub('\\sIndependent\\sCity', '\\sCity', county_ec$county)
county_ec$county <- gsub('\\sStaunton\\sWaynesboro', '', county_ec$county)

# manual fix for counties with area = 0
# scrape areas from Wikipedia, as Census data was inadequate
zeroarea <- county_ec[which(county_ec$area == 0), -c(2:7)] %>% 
  unique()
zeroarea$wiki_link <- gsub(' City and Borough', '', zeroarea$county)
zeroarea$wiki_link <- gsub(' Municipality', '', zeroarea$wiki_link)
zeroarea$wiki_link <- gsub('[A-Z]{2}', '', 
                           paste('https://en.wikipedia.org/wiki/', 
                                 gsub('\\s', '_', zeroarea$wiki_link), ',_', 
                                 gsub('AK', 'Alaska', zeroarea$state), 
                                 gsub('CO', 'Colorado', zeroarea$state), 
                                 sep = ''))
for(i in 1:nrow(zeroarea)) {
  zeroarea$area[i] <- gsub(',', '', read_html(zeroarea$wiki_link[i]) %>% 
                        html_nodes(xpath = '//th[contains(., "Area")]//parent::tr/following-sibling::tr[1]/td') %>% 
                        html_text() %>% 
                        str_extract('[\\d\\.,]+'))[1] %>% 
  as.numeric()
  county_ec$area[which(county_ec$county == zeroarea$county[i])] <- zeroarea$area[which(zeroarea$county == zeroarea$county[i])]
}

county_data <- merge(cp, county_ec, 
                     by = c('county', 'state', 'year')) %>% 
  mutate(empl_rate = empl_count / total_count, 
         pop_dens = total_count / area)
# 'county_data' to wide
st <- list()
for (i in 2013:2016) {
  df <- county_data %>% 
    filter(year == i)
  df <- df[, -3]
  df_names <- names(df)
  df_names[3] <- paste('p', i %>% as.character(), sep = '')
  df_names[4] <- paste('pop', i %>% as.character(), sep = '')
  df_names[5] <- paste('pr', i %>% as.character(), sep = '')
  df_names[6] <- paste('ec', i %>% as.character(), sep = '')
  df_names[7] <- paste('f_inc', i %>% as.character(), sep = '')
  df_names[8] <- paste('inc', i %>% as.character(), sep = '')
  df_names[9] <- paste('nf_inc', i %>% as.character(), sep = '')
  df_names[10] <- paste('pf_inc', i %>% as.character(), sep = '')
  df_names[12] <- paste('er', i %>% as.character(), sep = '')
  df_names[13] <- paste('pop_dens', i %>% as.character(), sep = '')
  names(df) <- df_names
  st[[i - 2012]] <- df
}
county_data_wide <- merge(st[[1]], st[[2]], 
                          by = c('state', 'county', 'area'))
county_data_wide <- merge(county_data_wide, st[[3]], 
                          by = c('state', 'county', 'area'))
county_data_wide <- merge(county_data_wide, st[[4]], 
                          by = c('state', 'county', 'area'))

# data set of poverty, by state
state_poverty <- table(county_data_wide$state) %>% 
  data.frame()
names(state_poverty) <- c('state', 'n_counties')
# 51 = number of states, including Washington D.C.
# create a data table of summary stats per state
sp <- list()
for(i in 1:51) {
  state_stats <- data.frame(summary(county_data_wide[which(county_data_wide$state == state_poverty$state[i]), ]))[13:258, 2:3] %>% 
  mutate(stat = gsub('\\s*\\:', '', 
                     str_extract(Freq, '.*\\:')), 
         Freq = gsub(':', '', 
                     str_extract(Freq, '\\:.+')) %>% 
           as.numeric()) %>% 
  reshape(idvar = 'Var2', timevar = 'stat', direction = 'wide') %>%
  mutate(Year = str_extract(Var2, '\\d{4}'), 
         Var2 = str_extract(Var2, '[a-z_]{1,7}'), 
         state = state_poverty$state[i])
  sp[[i]] <- state_stats
}
# clean data, rename variables and reorganize column order
state_poverty <- merge(state_poverty, do.call(rbind, sp), by = 'state') %>% 
  select(sort(current_vars()))
state_poverty <- state_poverty[, c(8, 10, 7, 9, 1:6)]
state_poverty$Var2[state_poverty$Var2 == 'p'] <- 'poverty_count'
state_poverty$Var2[state_poverty$Var2 == 'pr'] <- 'poverty_rate'
state_poverty$Var2[state_poverty$Var2 == 'pop'] <- 'total_count'
state_poverty$Var2[state_poverty$Var2 == 'ec'] <- 'empl_count'
state_poverty$Var2[state_poverty$Var2 == 'f_inc'] <- 'farm_income'
state_poverty$Var2[state_poverty$Var2 == 'inc'] <- 'income'
state_poverty$Var2[state_poverty$Var2 == 'f_inc'] <- 'farm_income'
state_poverty$Var2[state_poverty$Var2 == 'nf_inc'] <- 'nonfarm_income'
state_poverty$Var2[state_poverty$Var2 == 'pf_inc'] <- 'p_farm_income'
state_poverty$Var2[state_poverty$Var2 == 'er'] <- 'empl_rate'
state_poverty$Var2[state_poverty$Var2 == 'area'] <- 'area'
state_poverty$Var2[state_poverty$Var2 == 'pop_den'] <- 'pop_dens'
names(state_poverty) <- c('state', 'year', 'n_counties', 'countystat', '1st_quantile', '3rd_quantile', 'max', 'mean', 'median', 'minimum')
# replace NAs for 'year'
for(i in 1:nrow(state_poverty)) {
  if(state_poverty$year[i] %>% is.na) {
    state_poverty$year[i] <- state_poverty$year[i + 1]
  }
}
# create variables measuring the standard deviation and variance of county stats, per year
state_variance <- c()
state_sd <- c()
for(i in 1:nrow(state_poverty)) {
  state_variance[i] <- county_data %>% 
    filter(state == state_poverty$state[i] & year == state_poverty$year[i]) %>% 
    select(state_poverty$countystat[i]) %>% 
    var()
  state_sd[i] <- sapply(county_data %>% 
                          filter(state == state_poverty$state[i] & year == state_poverty$year[i]) %>% 
                          select(state_poverty$countystat[i]), 
                        sd)[[1]][1]
}
state_poverty$var <- state_variance
state_poverty$sd <- state_sd
state_poverty[is.na(state_poverty)] <- 0
# merge state stats with statistics on total state population ('state_ec'), total state poverty count, and overall state poverty rate (aggregating 'cp')
# change 'state_ec' names from full words to 2-letter abbreviations
for(i in 1:nrow(state_ec)) {
  if(state_ec$state[i] == 'District of Columbia') {
    state_ec$state[i] <- 'DC'
  } else {
    state_ec$state[i] <- state.abb[grep(state_ec$state[i], state.name)]
  }
}

number of items to replace is not a multiple of replacement lengthnumber of items to replace is not a multiple of replacement lengthnumber of items to replace is not a multiple of replacement lengthnumber of items to replace is not a multiple of replacement length

state_ec <- merge(state_ec, 
                  merge(aggregate(poverty_count ~ state + year, cp, sum), 
                        aggregate(total_count ~ state + year, cp, sum), 
                        by = c('state', 'year')), 
                  by = c('state', 'year')) %>% 
  mutate(poverty_rate = poverty_count / total_count, 
         empl_rate = empl_count / total_count)
state_ec <- merge(state_poverty, state_ec, by = c('state', 'year'))
# rename state-wide variables for clarity
names(state_ec) <- c('state', 'year', 'n_counties', 
                     'countystat', '1st_quantile', '3rd_quantile', 
                     'max', 'mean', 'median', 
                     'minimum', 'var', 'sd', 
                     'state_empl_count', 'state_farm_income', 'state_avg_income', 
                     'state_nonfarm_income', 'state_p_farm_income', 'state_area', 
                     'state_poverty_count', 'state_total_count', 'state_poverty_rate', 
                     'state_empl_rate')

# per-state crime data, merge with 'state_data'
state_crime <- read.csv(file = 'estimated_crimes.csv', header = TRUE, sep = ',')
state_data <- merge(state_ec, state_crime[, -c(2, 4:5, 15)], 
                    by.x = c('year', 'state'), by.y = c('year', 'state_abbr')) %>% 
  mutate(state_vcrime_rate = violent_crime / state_total_count,
         state_pcrime_rate = property_crime / state_total_count,
         state_pop_dens = state_total_count / state_area)
state_data <- state_data[, c(1:3, 20, 18, 34, 15, 4:12, 17, 13, 22, 19, 21, 23, 32, 28, 33)]
state_data$year <- state_data$year %>% as.numeric()
# explain 'state_data' variables below
# year: 2013-2016
# state: state, within the United States
# n_counties: number of counties in this state
# state_total_count: state population
# countystat: type of statistic, measured within counties of each state
## empl_count: number of employed individuals working per county
## farm_income: amount of income from farming per county
## nonfarm_income: amount of income from non-farming industries per county
## p_farm_income: ratio of 'farm_income' to 'nonfarm_income' per county, can be negative, indicator of how rural counties are
## empl_rate: ratio of employed individuals working per county, can be >1
## poverty_rate: percentage of individuals in poverty per county
## poverty_count: number of individuals in poverty per county
## total_count: population per county
# 1st_quantile:sd: statistics on counties, per state and year
# state_avg_income: average income per individual per state and year
# state_p_farm_income: ratio of 'farm_income' to 'nonfarm_income' per state and year, can be negative, indicator of how rural counties are
# state_poverty_count: number of individuals in poverty per state and year
# state_poverty_rate: percentage of individuals in poverty per state and year
# state_empl_count: number of employed individuals working per state and year
# state_empl_rate: ratio of employed individuals working per state and year
# violent_crime: number of violent crimes per state and year
# state_vcrime_rate: rate of violent crimes per state and year
# property_crime: number of property crimes per state and year
# state_pcrime_rate: rate of property crimes per state and year
# save 'state_data' as .csv file
write.csv(state_data, file = 'state_data.csv')

# various regressions
lm(state_vcrime_rate ~ state + year + var, 
   data = state_data %>% filter(countystat == 'poverty_rate')) %>% 
  msummary()

               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept) -0.16600898  0.02987540  -5.557 0.000000121172815565 ***
stateAL     -0.00256400  0.00016705 -15.349 < 0.0000000000000002 ***
stateAR     -0.00185904  0.00020580  -9.033 0.000000000000000738 ***
stateAZ     -0.00282847  0.00017330 -16.321 < 0.0000000000000002 ***
stateCA     -0.00280151  0.00019896 -14.081 < 0.0000000000000002 ***
stateCO     -0.00389012  0.00017447 -22.297 < 0.0000000000000002 ***
stateCT     -0.00434224  0.00032718 -13.272 < 0.0000000000000002 ***
stateDC      0.00709309  0.00037155  19.090 < 0.0000000000000002 ***
stateDE     -0.00146553  0.00036157  -4.053 0.000080609750724398 ***
stateFL     -0.00233992  0.00020156 -11.609 < 0.0000000000000002 ***
stateGA     -0.00348890  0.00018308 -19.056 < 0.0000000000000002 ***
stateHI     -0.00444839  0.00020108 -22.123 < 0.0000000000000002 ***
stateIA     -0.00396951  0.00029232 -13.579 < 0.0000000000000002 ***
stateID     -0.00478225  0.00022306 -21.440 < 0.0000000000000002 ***
stateIL     -0.00296354  0.00022463 -13.193 < 0.0000000000000002 ***
stateIN     -0.00304924  0.00025604 -11.909 < 0.0000000000000002 ***
stateKS     -0.00313147  0.00024636 -12.711 < 0.0000000000000002 ***
stateKY     -0.00520906  0.00018352 -28.384 < 0.0000000000000002 ***
stateLA     -0.00191695  0.00017918 -10.699 < 0.0000000000000002 ***
stateMA     -0.00281371  0.00026630 -10.566 < 0.0000000000000002 ***
stateMD     -0.00233531  0.00018345 -12.730 < 0.0000000000000002 ***
stateME     -0.00553603  0.00030059 -18.417 < 0.0000000000000002 ***
stateMI     -0.00248043  0.00024486 -10.130 < 0.0000000000000002 ***
stateMN     -0.00444102  0.00028297 -15.694 < 0.0000000000000002 ***
stateMO     -0.00219310  0.00020779 -10.555 < 0.0000000000000002 ***
stateMS     -0.00474653  0.00021443 -22.135 < 0.0000000000000002 ***
stateMT     -0.00375591  0.00017411 -21.572 < 0.0000000000000002 ***
stateNC     -0.00346685  0.00020981 -16.524 < 0.0000000000000002 ***
stateND     -0.00459135  0.00017237 -26.637 < 0.0000000000000002 ***
stateNE     -0.00403924  0.00026724 -15.115 < 0.0000000000000002 ***
stateNH     -0.00463935  0.00033665 -13.781 < 0.0000000000000002 ***
stateNJ     -0.00434025  0.00023178 -18.725 < 0.0000000000000002 ***
stateNM     -0.00078217  0.00016708  -4.681 0.000006292790458193 ***
stateNV     -0.00010532  0.00026927  -0.391              0.69625    
stateNY     -0.00301692  0.00025059 -12.039 < 0.0000000000000002 ***
stateOH     -0.00400565  0.00022659 -17.678 < 0.0000000000000002 ***
stateOK     -0.00250750  0.00022000 -11.398 < 0.0000000000000002 ***
stateOR     -0.00418656  0.00028998 -14.437 < 0.0000000000000002 ***
statePA     -0.00356266  0.00029106 -12.240 < 0.0000000000000002 ***
stateRI     -0.00448894  0.00024850 -18.064 < 0.0000000000000002 ***
stateSC     -0.00183783  0.00019820  -9.273 < 0.0000000000000002 ***
stateSD     -0.00445006  0.00053285  -8.351 0.000000000000040245 ***
stateTN     -0.00062988  0.00023868  -2.639              0.00919 ** 
stateTX     -0.00292348  0.00016753 -17.450 < 0.0000000000000002 ***
stateUT     -0.00470278  0.00019850 -23.692 < 0.0000000000000002 ***
stateVA     -0.00509680  0.00017780 -28.666 < 0.0000000000000002 ***
stateVT     -0.00553576  0.00032929 -16.811 < 0.0000000000000002 ***
stateWA     -0.00399023  0.00023101 -17.273 < 0.0000000000000002 ***
stateWI     -0.00389852  0.00026240 -14.857 < 0.0000000000000002 ***
stateWV     -0.00362839  0.00022946 -15.812 < 0.0000000000000002 ***
stateWY     -0.00481323  0.00021980 -21.898 < 0.0000000000000002 ***
year         0.00008573  0.00001484   5.775 0.000000042337838656 ***
var          0.17040076  0.08546137   1.994              0.04796 *  

Residual standard error: 0.0002361 on 151 degrees of freedom
Multiple R-squared:  0.9892,    Adjusted R-squared:  0.9855 
F-statistic: 267.2 on 52 and 151 DF,  p-value: < 0.00000000000000022

lm(state_vcrime_rate ~ state + year + var, 
   data = state_data %>% filter(countystat == 'pop_dens')) %>% 
  msummary()

                    Estimate       Std. Error t value             Pr(>|t|)    
(Intercept) -0.1827872107075  0.0303324801240  -6.026     0.00000001231859 ***
stateAL     -0.0025493182041  0.0001668002380 -15.284 < 0.0000000000000002 ***
stateAR     -0.0020979024549  0.0001667972137 -12.578 < 0.0000000000000002 ***
stateAZ     -0.0029184502459  0.0001668019413 -17.497 < 0.0000000000000002 ***
stateCA     -0.0028741075426  0.0001807304400 -15.903 < 0.0000000000000002 ***
stateCO     -0.0039084943962  0.0001715695923 -22.781 < 0.0000000000000002 ***
stateCT     -0.0048580721635  0.0001682370770 -28.876 < 0.0000000000000002 ***
stateDC      0.0064310304350  0.0001667963389  38.556 < 0.0000000000000002 ***
stateDE     -0.0020356735023  0.0001701650369 -11.963 < 0.0000000000000002 ***
stateFL     -0.0025321098431  0.0001675660907 -15.111 < 0.0000000000000002 ***
stateGA     -0.0033014425537  0.0001677913113 -19.676 < 0.0000000000000002 ***
stateHI     -0.0046580265936  0.0001669320526 -27.904 < 0.0000000000000002 ***
stateIA     -0.0044460794373  0.0001667988342 -26.655 < 0.0000000000000002 ***
stateID     -0.0050759510004  0.0001667974490 -30.432 < 0.0000000000000002 ***
stateIL     -0.0032158361676  0.0001683788492 -19.099 < 0.0000000000000002 ***
stateIN     -0.0034169488730  0.0001670618208 -20.453 < 0.0000000000000002 ***
stateKS     -0.0034863016927  0.0001668255622 -20.898 < 0.0000000000000002 ***
stateKY     -0.0050451230154  0.0001668970430 -30.229 < 0.0000000000000002 ***
stateLA     -0.0017758906576  0.0001668862211 -10.641 < 0.0000000000000002 ***
stateMA     -0.0026129656517  0.0003419663795  -7.641     0.00000000000231 ***
stateMD     -0.0020062377339  0.0002870198856  -6.990     0.00000000008269 ***
stateME     -0.0060332541499  0.0001667973205 -36.171 < 0.0000000000000002 ***
stateMI     -0.0028024760819  0.0001676667651 -16.715 < 0.0000000000000002 ***
stateMN     -0.0048562286115  0.0001679446942 -28.916 < 0.0000000000000002 ***
stateMO     -0.0023798550820  0.0001693199909 -14.055 < 0.0000000000000002 ***
stateMS     -0.0044775551126  0.0001667965858 -26.844 < 0.0000000000000002 ***
stateMT     -0.0038547370663  0.0001667963284 -23.110 < 0.0000000000000002 ***
stateNC     -0.0037039233242  0.0001669848897 -22.181 < 0.0000000000000002 ***
stateND     -0.0045058887798  0.0001667963239 -27.014 < 0.0000000000000002 ***
stateNE     -0.0044474679715  0.0001668400463 -26.657 < 0.0000000000000002 ***
stateNH     -0.0052165157168  0.0001668195873 -31.270 < 0.0000000000000002 ***
stateNJ     -0.0029459452911  0.0008497418100  -3.467             0.000686 ***
stateNM     -0.0007657653601  0.0001668009087  -4.591     0.00000922818479 ***
stateNV     -0.0005257684625  0.0001667967761  -3.152             0.001955 ** 
stateNY      0.0112834420635  0.0071310049495   1.582             0.115672    
stateOH     -0.0042770322678  0.0001676181947 -25.517 < 0.0000000000000002 ***
stateOK     -0.0027859216144  0.0001668337819 -16.699 < 0.0000000000000002 ***
stateOR     -0.0046375660160  0.0001671314123 -27.748 < 0.0000000000000002 ***
statePA     -0.0035689816481  0.0002824186807 -12.637 < 0.0000000000000002 ***
stateRI     -0.0047882652724  0.0001700080362 -28.165 < 0.0000000000000002 ***
stateSC     -0.0020469628469  0.0001668071027 -12.271 < 0.0000000000000002 ***
stateSD     -0.0034410995171  0.0001667963070 -20.631 < 0.0000000000000002 ***
stateTN     -0.0009597260154  0.0001668712505  -5.751     0.00000004759584 ***
stateTX     -0.0029277776899  0.0001671974513 -17.511 < 0.0000000000000002 ***
stateUT     -0.0048991442243  0.0001670196539 -29.333 < 0.0000000000000002 ***
stateVA     -0.0046485052382  0.0003234171580 -14.373 < 0.0000000000000002 ***
stateVT     -0.0061011334634  0.0001667965308 -36.578 < 0.0000000000000002 ***
stateWA     -0.0042988167742  0.0001668642171 -25.762 < 0.0000000000000002 ***
stateWI     -0.0042972375386  0.0001668136639 -25.761 < 0.0000000000000002 ***
stateWV     -0.0039400637857  0.0001667998460 -23.622 < 0.0000000000000002 ***
stateWY     -0.0050985128235  0.0001667963351 -30.567 < 0.0000000000000002 ***
year         0.0000943869796  0.0000150569724   6.269     0.00000000362896 ***
var         -0.0000000002637  0.0000000001281  -2.058             0.041290 *  

Residual standard error: 0.0002359 on 151 degrees of freedom
Multiple R-squared:  0.9893,    Adjusted R-squared:  0.9856 
F-statistic: 267.6 on 52 and 151 DF,  p-value: < 0.00000000000000022

lm(state_vcrime_rate ~ state + year + mean, 
   data = state_data %>% filter(countystat == 'p_farm_income')) %>% 
  msummary()

               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept) -0.12103252  0.03431532  -3.527             0.000557 ***
stateAL     -0.00244393  0.00016927 -14.438 < 0.0000000000000002 ***
stateAR     -0.00197214  0.00017094 -11.537 < 0.0000000000000002 ***
stateAZ     -0.00289953  0.00016505 -17.568 < 0.0000000000000002 ***
stateCA     -0.00294148  0.00016706 -17.607 < 0.0000000000000002 ***
stateCO     -0.00393034  0.00016628 -23.636 < 0.0000000000000002 ***
stateCT     -0.00490279  0.00016487 -29.738 < 0.0000000000000002 ***
stateDC      0.00643088  0.00016487  39.007 < 0.0000000000000002 ***
stateDE     -0.00206329  0.00016553 -12.465 < 0.0000000000000002 ***
stateFL     -0.00249554  0.00016672 -14.969 < 0.0000000000000002 ***
stateGA     -0.00324939  0.00016792 -19.351 < 0.0000000000000002 ***
stateHI     -0.00467270  0.00016487 -28.343 < 0.0000000000000002 ***
stateIA     -0.00425524  0.00017857 -23.830 < 0.0000000000000002 ***
stateID     -0.00487300  0.00018018 -27.046 < 0.0000000000000002 ***
stateIL     -0.00314827  0.00016987 -18.534 < 0.0000000000000002 ***
stateIN     -0.00335816  0.00016720 -20.085 < 0.0000000000000002 ***
stateKS     -0.00327491  0.00018218 -17.976 < 0.0000000000000002 ***
stateKY     -0.00500224  0.00016602 -30.131 < 0.0000000000000002 ***
stateLA     -0.00173208  0.00016603 -10.433 < 0.0000000000000002 ***
stateMA     -0.00322817  0.00016487 -19.581 < 0.0000000000000002 ***
stateMD     -0.00245627  0.00016523 -14.866 < 0.0000000000000002 ***
stateME     -0.00602590  0.00016489 -36.544 < 0.0000000000000002 ***
stateMI     -0.00282840  0.00016490 -17.152 < 0.0000000000000002 ***
stateMN     -0.00477015  0.00017090 -27.913 < 0.0000000000000002 ***
stateMO     -0.00236022  0.00016728 -14.109 < 0.0000000000000002 ***
stateMS     -0.00443108  0.00016572 -26.739 < 0.0000000000000002 ***
stateMT     -0.00367647  0.00017663 -20.814 < 0.0000000000000002 ***
stateNC     -0.00364337  0.00016712 -21.801 < 0.0000000000000002 ***
stateND     -0.00436628  0.00017217 -25.360 < 0.0000000000000002 ***
stateNE     -0.00386418  0.00026730 -14.456 < 0.0000000000000002 ***
stateNH     -0.00522206  0.00016487 -31.675 < 0.0000000000000002 ***
stateNJ     -0.00465974  0.00016487 -28.264 < 0.0000000000000002 ***
stateNM     -0.00066467  0.00016894  -3.934             0.000127 ***
stateNV     -0.00055989  0.00016529  -3.387             0.000900 ***
stateNY     -0.00336834  0.00016504 -20.409 < 0.0000000000000002 ***
stateOH     -0.00429453  0.00016497 -26.032 < 0.0000000000000002 ***
stateOK     -0.00266827  0.00017076 -15.626 < 0.0000000000000002 ***
stateOR     -0.00457832  0.00016737 -27.355 < 0.0000000000000002 ***
statePA     -0.00402305  0.00016495 -24.389 < 0.0000000000000002 ***
stateRI     -0.00485617  0.00016487 -29.455 < 0.0000000000000002 ***
stateSC     -0.00203865  0.00016492 -12.361 < 0.0000000000000002 ***
stateSD     -0.00303679  0.00021883 -13.877 < 0.0000000000000002 ***
stateTN     -0.00096544  0.00016487  -5.856         0.0000000286 ***
stateTX     -0.00284826  0.00016892 -16.862 < 0.0000000000000002 ***
stateUT     -0.00480114  0.00016994 -28.253 < 0.0000000000000002 ***
stateVA     -0.00520559  0.00016493 -31.562 < 0.0000000000000002 ***
stateVT     -0.00607852  0.00016507 -36.824 < 0.0000000000000002 ***
stateWA     -0.00424245  0.00016654 -25.474 < 0.0000000000000002 ***
stateWI     -0.00425670  0.00016566 -25.696 < 0.0000000000000002 ***
stateWV     -0.00394624  0.00016487 -23.935 < 0.0000000000000002 ***
stateWY     -0.00507286  0.00016511 -30.723 < 0.0000000000000002 ***
year         0.00006373  0.00001703   3.741             0.000259 ***
mean        -0.00020643  0.00007347  -2.810             0.005616 ** 

Residual standard error: 0.0002332 on 151 degrees of freedom
Multiple R-squared:  0.9895,    Adjusted R-squared:  0.9859 
F-statistic:   274 on 52 and 151 DF,  p-value: < 0.00000000000000022

lm(state_vcrime_rate ~ state + year + var, 
   data = state_data %>% filter(countystat == 'empl_rate')) %>% 
  msummary()

               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept) -0.16584128  0.03046629  -5.443  0.00000020690236819 ***
stateAL     -0.00232649  0.00028211  -8.247  0.00000000000007374 ***
stateAR     -0.00187305  0.00028276  -6.624  0.00000000057725935 ***
stateAZ     -0.00269917  0.00027963  -9.653 < 0.0000000000000002 ***
stateCA     -0.00281504  0.00026403 -10.662 < 0.0000000000000002 ***
stateCO     -0.00383766  0.00022849 -16.796 < 0.0000000000000002 ***
stateCT     -0.00467380  0.00028559 -16.365 < 0.0000000000000002 ***
stateDC      0.00668593  0.00030639  21.822 < 0.0000000000000002 ***
stateDE     -0.00187060  0.00028960  -6.459  0.00000000136363296 ***
stateFL     -0.00235305  0.00027167  -8.662  0.00000000000000661 ***
stateGA     -0.00318462  0.00022905 -13.903 < 0.0000000000000002 ***
stateHI     -0.00442819  0.00029723 -14.898 < 0.0000000000000002 ***
stateIA     -0.00422026  0.00028415 -14.852 < 0.0000000000000002 ***
stateID     -0.00520768  0.00021352 -24.390 < 0.0000000000000002 ***
stateIL     -0.00304324  0.00027795 -10.949 < 0.0000000000000002 ***
stateIN     -0.00322055  0.00027460 -11.728 < 0.0000000000000002 ***
stateKS     -0.00328264  0.00027010 -12.153 < 0.0000000000000002 ***
stateKY     -0.00485795  0.00026157 -18.572 < 0.0000000000000002 ***
stateLA     -0.00160218  0.00025089  -6.386  0.00000000199158404 ***
stateMA     -0.00307112  0.00023035 -13.332 < 0.0000000000000002 ***
stateMD     -0.00226145  0.00028240  -8.008  0.00000000000028970 ***
stateME     -0.00580699  0.00028394 -20.451 < 0.0000000000000002 ***
stateMI     -0.00261719  0.00027827  -9.405 < 0.0000000000000002 ***
stateMN     -0.00467744  0.00027726 -16.870 < 0.0000000000000002 ***
stateMO     -0.00222256  0.00027575  -8.060  0.00000000000021527 ***
stateMS     -0.00427042  0.00026828 -15.918 < 0.0000000000000002 ***
stateMT     -0.00364065  0.00027312 -13.330 < 0.0000000000000002 ***
stateNC     -0.00351146  0.00026909 -13.049 < 0.0000000000000002 ***
stateND     -0.00435704  0.00022522 -19.345 < 0.0000000000000002 ***
stateNE     -0.00422841  0.00028351 -14.914 < 0.0000000000000002 ***
stateNH     -0.00499346  0.00028503 -17.519 < 0.0000000000000002 ***
stateNJ     -0.00443768  0.00028050 -15.820 < 0.0000000000000002 ***
stateNM     -0.00056434  0.00026533  -2.127              0.03505 *  
stateNV     -0.00122579  0.00072234  -1.697              0.09176 .  
stateNY     -0.00323640  0.00022823 -14.180 < 0.0000000000000002 ***
stateOH     -0.00409035  0.00027859 -14.682 < 0.0000000000000002 ***
stateOK     -0.00258715  0.00026694  -9.692 < 0.0000000000000002 ***
stateOR     -0.00443097  0.00028468 -15.565 < 0.0000000000000002 ***
statePA     -0.00382159  0.00027513 -13.890 < 0.0000000000000002 ***
stateRI     -0.00461594  0.00029420 -15.690 < 0.0000000000000002 ***
stateSC     -0.00183411  0.00027537  -6.660  0.00000000047726822 ***
stateSD     -0.00324702  0.00025773 -12.598 < 0.0000000000000002 ***
stateTN     -0.00075813  0.00027152  -2.792              0.00591 ** 
stateTX     -0.00277209  0.00024685 -11.230 < 0.0000000000000002 ***
stateUT     -0.00472267  0.00025786 -18.315 < 0.0000000000000002 ***
stateVA     -0.00503090  0.00025307 -19.879 < 0.0000000000000002 ***
stateVT     -0.00590757  0.00025777 -22.918 < 0.0000000000000002 ***
stateWA     -0.00408172  0.00028349 -14.398 < 0.0000000000000002 ***
stateWI     -0.00407603  0.00028290 -14.408 < 0.0000000000000002 ***
stateWV     -0.00373564  0.00026742 -13.969 < 0.0000000000000002 ***
stateWY     -0.00490828  0.00025467 -19.273 < 0.0000000000000002 ***
year         0.00008585  0.00001514   5.669  0.00000007092442572 ***
var          0.00243152  0.00244261   0.995              0.32111    

Residual standard error: 0.0002384 on 151 degrees of freedom
Multiple R-squared:  0.989, Adjusted R-squared:  0.9853 
F-statistic:   262 on 52 and 151 DF,  p-value: < 0.00000000000000022

lm(state_vcrime_rate ~ state + year + mean, 
   data = state_data %>% filter(countystat == 'income')) %>% 
    msummary()

                  Estimate     Std. Error t value             Pr(>|t|)    
(Intercept) -0.30803238818  0.04063145024  -7.581     0.00000000000323 ***
stateAL     -0.00416643472  0.00037861672 -11.004 < 0.0000000000000002 ***
stateAR     -0.00379390332  0.00039419386  -9.624 < 0.0000000000000002 ***
stateAZ     -0.00455643977  0.00038257297 -11.910 < 0.0000000000000002 ***
stateCA     -0.00361869746  0.00020342279 -17.789 < 0.0000000000000002 ***
stateCO     -0.00480949741  0.00023529420 -20.440 < 0.0000000000000002 ***
stateCT     -0.00477181893  0.00016045448 -29.739 < 0.0000000000000002 ***
stateDC      0.00747919334  0.00027353978  27.342 < 0.0000000000000002 ***
stateDE     -0.00300575797  0.00024859396 -12.091 < 0.0000000000000002 ***
stateFL     -0.00391377937  0.00032794047 -11.934 < 0.0000000000000002 ***
stateGA     -0.00502341846  0.00039215544 -12.810 < 0.0000000000000002 ***
stateHI     -0.00562997717  0.00025814702 -21.809 < 0.0000000000000002 ***
stateIA     -0.00529545393  0.00023994099 -22.070 < 0.0000000000000002 ***
stateID     -0.00642844887  0.00032842199 -19.574 < 0.0000000000000002 ***
stateIL     -0.00439976980  0.00028915377 -15.216 < 0.0000000000000002 ***
stateIN     -0.00472026937  0.00031592791 -14.941 < 0.0000000000000002 ***
stateKS     -0.00435004386  0.00024152096 -18.011 < 0.0000000000000002 ***
stateKY     -0.00668803845  0.00038176090 -17.519 < 0.0000000000000002 ***
stateLA     -0.00313738740  0.00032823436  -9.558 < 0.0000000000000002 ***
stateMA     -0.00291909976  0.00017110362 -17.060 < 0.0000000000000002 ***
stateMD     -0.00296903649  0.00018844457 -15.755 < 0.0000000000000002 ***
stateME     -0.00721846172  0.00029767852 -24.249 < 0.0000000000000002 ***
stateMI     -0.00424793720  0.00033951924 -12.512 < 0.0000000000000002 ***
stateMN     -0.00575253856  0.00024130247 -23.840 < 0.0000000000000002 ***
stateMO     -0.00397575572  0.00036342436 -10.940 < 0.0000000000000002 ***
stateMS     -0.00619859958  0.00039918729 -15.528 < 0.0000000000000002 ***
stateMT     -0.00493053953  0.00027843616 -17.708 < 0.0000000000000002 ***
stateNC     -0.00519079872  0.00035092482 -14.792 < 0.0000000000000002 ***
stateND     -0.00483040704  0.00017247324 -28.007 < 0.0000000000000002 ***
stateNE     -0.00492962454  0.00018754939 -26.284 < 0.0000000000000002 ***
stateNH     -0.00568886127  0.00018667213 -30.475 < 0.0000000000000002 ***
stateNJ     -0.00469149637  0.00015812475 -29.670 < 0.0000000000000002 ***
stateNM     -0.00217797172  0.00033938978  -6.417     0.00000000169371 ***
stateNV     -0.00164746026  0.00028636598  -5.753     0.00000004720661 ***
stateNY     -0.00426064475  0.00024374648 -17.480 < 0.0000000000000002 ***
stateOH     -0.00557217730  0.00031170825 -17.876 < 0.0000000000000002 ***
stateOK     -0.00410503339  0.00032108840 -12.785 < 0.0000000000000002 ***
stateOR     -0.00590509273  0.00030890806 -19.116 < 0.0000000000000002 ***
statePA     -0.00508819295  0.00027392141 -18.575 < 0.0000000000000002 ***
stateRI     -0.00494336731  0.00015908400 -31.074 < 0.0000000000000002 ***
stateSC     -0.00365682220  0.00037691226  -9.702 < 0.0000000000000002 ***
stateSD     -0.00420312439  0.00022655218 -18.553 < 0.0000000000000002 ***
stateTN     -0.00251144894  0.00036447525  -6.891     0.00000000014088 ***
stateTX     -0.00404450340  0.00028141340 -14.372 < 0.0000000000000002 ***
stateUT     -0.00633218074  0.00034044793 -18.600 < 0.0000000000000002 ***
stateVA     -0.00626660696  0.00027351357 -22.912 < 0.0000000000000002 ***
stateVT     -0.00689741504  0.00023175313 -29.762 < 0.0000000000000002 ***
stateWA     -0.00535514155  0.00027329230 -19.595 < 0.0000000000000002 ***
stateWI     -0.00533107190  0.00027023299 -19.728 < 0.0000000000000002 ***
stateWV     -0.00557482269  0.00038205970 -14.591 < 0.0000000000000002 ***
stateWY     -0.00537930721  0.00016895514 -31.839 < 0.0000000000000002 ***
year         0.00015847410  0.00002046534   7.744     0.00000000000130 ***
mean        -0.00000006697  0.00000001427  -4.693     0.00000598884176 ***

Residual standard error: 0.0002234 on 151 degrees of freedom
Multiple R-squared:  0.9904,    Adjusted R-squared:  0.9871 
F-statistic: 298.6 on 52 and 151 DF,  p-value: < 0.00000000000000022

state_data %>% 
  filter(countystat == 'poverty_rate') %>% 
  select(var) %>% 
  unique() %>% 
  ggplot(aes(x = var)) + 
  geom_density()

# decision trees
rpart(state_vcrime_rate ~ state + year, 
      data = state_data)

n= 2091 

node), split, n, deviance, yval
      * denotes terminal node

 1) root 2091 0.008022991000 0.003946431  
   2) state=AL,AR,AZ,CA,CO,CT,DE,FL,GA,HI,IA,ID,IL,IN,KS,KY,LA,MA,MD,ME,MI,MN,MO,MS,MT,NC,ND,NE,NH,NJ,NY,OH,OK,OR,PA,RI,SC,SD,TX,UT,VA,VT,WA,WI,WV,WY 1886 0.002262534000 0.003485893  
     4) state=CO,CT,HI,IA,ID,KY,ME,MN,MS,MT,ND,NE,NH,NJ,OH,OR,PA,RI,UT,VA,VT,WA,WI,WV,WY 1025 0.000337117100 0.002634185  
       8) state=CT,ID,KY,ME,MN,NH,RI,UT,VA,VT,WY 451 0.000085119900 0.002139351 *
       9) state=CO,HI,IA,MS,MT,ND,NE,NJ,OH,OR,PA,WA,WI,WV 574 0.000054796770 0.003022983 *
     5) state=AL,AR,AZ,CA,DE,FL,GA,IL,IN,KS,LA,MA,MD,MI,MO,NC,NY,OK,SC,SD,TX 861 0.000296705500 0.004499832  
      10) state=AZ,CA,GA,IL,IN,KS,MA,MI,NC,NY,OK,SD,TX 533 0.000064932130 0.004135954 *
      11) state=AL,AR,DE,FL,LA,MD,MO,SC 328 0.000046519510 0.005091133 *
   3) state=AK,DC,NM,NV,TN 205 0.001680338000 0.008183381  
     6) state=AK,NM,NV,TN 164 0.000055436340 0.006779060 *
     7) state=DC 41 0.000007764723 0.013800670 *

rpart(state_vcrime_rate ~ var, 
      data = state_data %>% filter(countystat == 'poverty_rate'))

n= 204 

node), split, n, deviance, yval
      * denotes terminal node

 1) root 204 0.000782781900 0.003948090  
   2) var>=0.0001590567 197 0.000380248900 0.003728483  
     4) var< 0.001100494 25 0.000020637570 0.002282980 *
     5) var>=0.001100494 172 0.000299781700 0.003938586  
      10) var< 0.003312984 126 0.000180747600 0.003761335 *
      11) var>=0.003312984 46 0.000104232200 0.004424099  
        22) var>=0.004113813 26 0.000039203400 0.003755888 *
        23) var< 0.004113813 20 0.000038327790 0.005292773  
          46) var< 0.003699606 8 0.000001079445 0.004317887 *
          47) var>=0.003699606 12 0.000024576320 0.005942696 *
   3) var< 0.0001590567 7 0.000125653500 0.010128460 *

rpart(state_vcrime_rate ~ ., 
      data = state_data)

n= 2091 

node), split, n, deviance, yval
      * denotes terminal node

 1) root 2091 0.008022991000 0.003946431  
   2) state=AL,AR,AZ,CA,CO,CT,DE,FL,GA,HI,IA,ID,IL,IN,KS,KY,LA,MA,MD,ME,MI,MN,MO,MS,MT,NC,ND,NE,NH,NJ,NY,OH,OK,OR,PA,RI,SC,SD,TX,UT,VA,VT,WA,WI,WV,WY 1886 0.002262534000 0.003485893  
     4) state=CO,CT,HI,IA,ID,KY,ME,MN,MS,MT,ND,NE,NH,NJ,OH,OR,PA,RI,UT,VA,VT,WA,WI,WV,WY 1025 0.000337117100 0.002634185  
       8) state=CT,ID,KY,ME,MN,NH,RI,UT,VA,VT,WY 451 0.000085119900 0.002139351 *
       9) state=CO,HI,IA,MS,MT,ND,NE,NJ,OH,OR,PA,WA,WI,WV 574 0.000054796770 0.003022983 *
     5) state=AL,AR,AZ,CA,DE,FL,GA,IL,IN,KS,LA,MA,MD,MI,MO,NC,NY,OK,SC,SD,TX 861 0.000296705500 0.004499832  
      10) state=AZ,CA,GA,IL,IN,KS,MA,MI,NC,NY,OK,SD,TX 533 0.000064932130 0.004135954 *
      11) state=AL,AR,DE,FL,LA,MD,MO,SC 328 0.000046519510 0.005091133 *
   3) state=AK,DC,NM,NV,TN 205 0.001680338000 0.008183381  
     6) state_avg_income< 62777.5 164 0.000055436340 0.006779060 *
     7) state_avg_income>=62777.5 41 0.000007764723 0.013800670 *

rpart(state_vcrime_rate ~ state + var, 
      data = state_data %>% filter(countystat == 'pop_dens'))

n= 204 

node), split, n, deviance, yval
      * denotes terminal node

 1) root 204 0.000782781900 0.003948090  
   2) state=AL,AR,AZ,CA,CO,CT,DE,FL,GA,HI,IA,ID,IL,IN,KS,KY,LA,MA,MD,ME,MI,MN,MO,MS,MT,NC,ND,NE,NH,NJ,NY,OH,OK,OR,PA,RI,SC,SD,TX,UT,VA,VT,WA,WI,WV,WY 184 0.000221201600 0.003487186  
     4) state=CO,CT,HI,IA,ID,KY,ME,MN,MS,MT,ND,NE,NH,NJ,OH,OR,PA,RI,UT,VA,VT,WA,WI,WV,WY 100 0.000033043870 0.002634751  
       8) state=CT,ID,KY,ME,MN,NH,RI,UT,VA,VT,WY 44 0.000008296347 0.002138552 *
       9) state=CO,HI,IA,MS,MT,ND,NE,NJ,OH,OR,PA,WA,WI,WV 56 0.000005402242 0.003024621 *
     5) state=AL,AR,AZ,CA,DE,FL,GA,IL,IN,KS,LA,MA,MD,MI,MO,NC,NY,OK,SC,SD,TX 84 0.000028987610 0.004501990  
      10) state=AZ,CA,GA,IL,IN,KS,MA,MI,NC,NY,OK,SD,TX 52 0.000006301939 0.004137407 *
      11) state=AL,AR,DE,FL,LA,MD,MO,SC 32 0.000004541995 0.005094437 *
   3) state=AK,DC,NM,NV,TN 20 0.000162886000 0.008188413  
     6) var>=811.4494 13 0.000004729091 0.006737831 *
     7) var< 811.4494 7 0.000080001360 0.010882350 *

plot(as.party(rpart(state_vcrime_rate ~ ., data = state_data)))

plot(as.party(rpart(state_vcrime_rate ~ var + state_avg_income, data = state_data %>% filter(countystat == 'income'))))

state_data_pca_df <- state_data[,c(1:7, 17:25)] %>% 
  unique()
rownames(state_data_pca_df) <- paste(state_data_pca_df$state, state_data_pca_df$year, sep = '')
state_data_pca <- state_data_pca_df[, -2] %>% 
  prcomp(scale = TRUE)
-1 * state_data_pca$rotation[, 1:4] %>% 
  round(3)

                       PC1    PC2    PC3    PC4
year                -0.002 -0.049  0.190 -0.593
n_counties          -0.246  0.111 -0.228  0.179
state_total_count   -0.419 -0.037  0.120  0.051
state_area          -0.088 -0.024 -0.051 -0.037
state_pop_dens       0.068 -0.529 -0.006  0.106
state_avg_income     0.037 -0.333  0.509  0.110
state_p_farm_income  0.076  0.141 -0.177  0.673
state_empl_count    -0.417 -0.047  0.138  0.062
state_empl_rate      0.131 -0.468  0.119  0.273
state_poverty_count -0.426 -0.039  0.050  0.034
state_poverty_rate  -0.137 -0.057 -0.550 -0.196
violent_crime       -0.419 -0.074  0.066  0.014
state_vcrime_rate   -0.018 -0.474 -0.245 -0.109
property_crime      -0.421 -0.049  0.022  0.047
state_pcrime_rate   -0.012 -0.338 -0.449 -0.066

summary(state_data_pca)

Importance of components:
                          PC1    PC2    PC3    PC4     PC5    PC6    PC7     PC8     PC9    PC10    PC11    PC12    PC13    PC14    PC15
Standard deviation     2.3236 1.7514 1.4704 1.0927 1.04345 0.9558 0.6842 0.55525 0.44079 0.34091 0.24546 0.13071 0.08364 0.05586 0.02395
Proportion of Variance 0.3599 0.2045 0.1442 0.0796 0.07259 0.0609 0.0312 0.02055 0.01295 0.00775 0.00402 0.00114 0.00047 0.00021 0.00004
Cumulative Proportion  0.3599 0.5644 0.7086 0.7882 0.86077 0.9217 0.9529 0.97343 0.98638 0.99413 0.99815 0.99929 0.99975 0.99996 1.00000

data.frame(sd = state_data_pca$sdev) %>% 
  rownames_to_column() %>% 
  mutate(rowname = parse_number(rowname), 
         totalVar = sum(state_data_pca$sdev^2), 
         pve = 100 * sd^2 / totalVar, 
         cusum = cumsum(pve)) %>% 
  ggplot(aes(x = rowname, y = pve)) + 
  geom_line()

data.frame(sd = state_data_pca$sdev) %>% 
  rownames_to_column() %>% 
  mutate(rowname = parse_number(rowname), 
         totalVar = sum(state_data_pca$sdev^2), 
         pve = 100 * sd^2 / totalVar, 
         cusum = cumsum(pve)) %>% 
  ggplot(aes(x = rowname, y = cusum)) + 
  geom_line()

pca_data <- state_data_pca$x %>% 
  data.frame() %>% 
  rownames_to_column() %>% 
  mutate(year = str_extract(rowname, '\\d+') %>% as.numeric(), 
         rowname = gsub('\\d{4}', '', rowname))
pca_data %>% 
  ggplot(aes(x = -PC1, y = -PC2, 
             colour = year %>% as.factor())) + 
  geom_text(aes(label = rowname), size = 3)

# averaged values, per year
aggregate(pca_data, by = list(pca_data$rowname), FUN = mean) %>% 
  ggplot(aes(x = -PC1, y = -PC2)) + 
  geom_text(aes(label = Group.1), size = 3)

argument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NA

# averaged values, per year
pca_data_alt <- aggregate(pca_data, by = list(pca_data$rowname), FUN = mean)

argument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NA

rownames(pca_data_alt) <- pca_data_alt$Group.1
pca_data_alt <- pca_data_alt[, 3:17]
pca_data_alt %>% 
  dist() %>% 
  hclust() %>% 
  plot(cex = 0.8)

pca_data2 <- state_data_pca$x %>% 
  data.frame() %>% 
  rownames_to_column() 
pca_data2 %>% 
  dist() %>% 
  hclust() %>% 
  plot(cex = 0.8, labels = pca_data2$rowname)

NAs introduced by coercion

# visualization
df <- aggregate(state_poverty_rate ~ state, data = state_data, mean)
df <- filter(df, 
             state != "DC") 
mycolors_True <- brewer.pal(9, "GnBu")
states <- tolower(rownames(USArrests))
map_50S <- map_data("state")
ggplot(data = df,
       aes(fill = state_poverty_rate)) +
  labs(x = "Longitude", y = "Latitude", title = "Average Poverty Rate in States") +
  geom_map(aes(map_id = states), map = map_50S) +
  expand_limits(x = map_50S$long, y = map_50S$lat) +
  theme_bw() + 
  scale_fill_gradientn(name = "Poverty Rate", 
                       colours = mycolors_True,
                       na.value = "grey") 

# change two-letter abbreviation to full state name for merging
abbr2state <- function(abbr){
  ab    <- tolower(c("AL",
             "AK", "AZ", "KS", "UT", "CO", "CT",
             "DE", "FL", "GA", "HI", "ID", "IL",
             "IN", "IA", "AR", "KY", "LA", "ME",
             "MD", "MA", "MI", "MN", "MS", "MO",
             "MT", "NE", "NV", "NH", "NJ", "NM",
             "NY", "NC", "ND", "OH", "OK", "OR",
             "PA", "RI", "SC", "SD", "TN", "TX",
             "CA", "VT", "VA", "WA", "WV", "WI",
             "WY", "DC"))
  st    <- c("Alabama",
             "Alaska", "Arizona", "Kansas",
             "Utah", "Colorado", "Connecticut",
             "Delaware", "Florida", "Georgia",
             "Hawaii", "Idaho", "Illinois",
             "Indiana", "Iowa", "Arkansas",
             "Kentucky", "Louisiana", "Maine",
             "Maryland", "Massachusetts", "Michigan",
             "Minnesota", "Mississippi", "Missouri",
             "Montana", "Nebraska", "Nevada",
             "New Hampshire", "New Jersey", "New Mexico",
             "New York", "North Carolina", "North Dakota",
             "Ohio", "Oklahoma", "Oregon",
             "Pennsylvania", "Rhode Island", "South Carolina",
             "South Dakota", "Tennessee", "Texas",
             "California", "Vermont", "Virginia",
             "Washington", "West Virginia", "Wisconsin",
             "Wyoming", "District of Columbia")
  st[match(tolower(abbr), ab)]
}
# visualization
# fixed map
merge(map_data("state"), 
                      aggregate(var ~ state, 
                                data = filter(state_data, countystat == "poverty_rate") %>% 
                                  select(year, state, var) %>% 
                                  unique(), 
                                FUN = mean) %>% 
                        mutate(state = abbr2state(state) %>% tolower()), 
                      by.x = 'region', by.y = 'state') %>% 
  ggplot() + 
  geom_polygon(aes(long, lat, 
                   group = group, 
                   fill = var)) + 
  scale_fill_gradientn(colours = mycolors_True, 
                       trans = 'log', 
                       name = 'Variance') + 
  theme_bw() + 
  labs(x = "Longitude", 
       y = "Latitude", 
       title = "Variance of Poverty Rate in Counties, Per State")

# visualization
merge(map_data("state"), 
                      aggregate(var ~ state, 
                                data = filter(state_data, countystat == "pop_dens") %>% 
                                  select(year, state, var) %>% 
                                  unique(), 
                                FUN = mean) %>% 
                        mutate(state = abbr2state(state) %>% tolower()), 
                      by.x = 'region', by.y = 'state') %>% 
  ggplot() + 
  geom_polygon(aes(long, lat, 
                   group = group, 
                   fill = var)) + 
  scale_fill_gradientn(colours = mycolors_True, 
                       trans = 'log', 
                       name = 'Variance') + 
  theme_bw() + 
  labs(x = "Longitude", 
       y = "Latitude", 
       title = "Variance of Population Density in Counties, Per State")

aggregate(state_pcrime_rate ~ state, data = state_data, mean) %>% 
  ggplot(aes(x = state, 
             y = state_pcrime_rate, 
             colour = state, 
             group = state)) +
  geom_label(aes(label = state), check_overlap = TRUE)+
  labs(x = "State", 
       y = "Average State Property Crime Rate Per State (including DC)", 
       title = "Average State Property Crime Rate Per State (including DC)")

Ignoring unknown parameters: check_overlap

df_Farm_Counties <- merge(aggregate(state_p_farm_income ~ state, data = state_data, mean), 
                          aggregate(n_counties ~ state, data = state_data, mean), 
                          by = c('state'))
rownames(df_Farm_Counties) <- df_Farm_Counties$state
set.seed(100)
fit <- kmeans(df_Farm_Counties[, 2:3], 5)
clusplot(df_Farm_Counties, 
         fit$cluster, 
         color = TRUE, 
         shade = TRUE, 
         labels = 2, 
         lines = 0,
         cex = 0.75,
         main = 'K-Means Clustering, Proportion of Farm Income \nand Number of Counties per State')

df_Farm_Poverty <- merge(aggregate(state_p_farm_income ~ state, data = state_data, mean), 
                         aggregate(state_poverty_rate ~ state, data = state_data, mean), 
                         by = c('state'))
rownames(df_Farm_Poverty) <- df_Farm_Poverty$state
set.seed(100)
fit2 <- kmeans(df_Farm_Poverty[, 2:3], 4)
clusplot(df_Farm_Poverty, 
         fit2$cluster, 
         color = TRUE, 
         shade = TRUE, 
         labels = 2, 
         lines = 0,
         cex = 0.75,
         main = 'K-Means Clustering, Proportion of Farm Income and Poverty Rate')

df_Crime_Poverty <- merge(aggregate(state_vcrime_rate ~ state, data = state_data, mean), 
                         aggregate(state_poverty_rate ~ state, data = state_data, mean), 
                         by = c('state'))
rownames(df_Crime_Poverty) <- df_Crime_Poverty$state
set.seed(100)
fit2 <- kmeans(df_Crime_Poverty[, 2:3], 4)
clusplot(df_Crime_Poverty, 
         fit2$cluster, 
         color = TRUE, 
         shade = TRUE, 
         labels = 2, 
         lines = 0,
         cex = 0.75,
         main = 'K-Means Clustering, Proportion of Violent Crime Rate and Poverty Rate')

stopCluster(myCluster)